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A NUMERICAL  ALGORITHM  FOR  REMOTE  SENSING  OF 


OCEAN  DENSITY  PROFILES  BY  ACOUSTIC  PULSES 

Y.  M.  Chen  and  D.  S.  Tsien 
Numerical  Computation  Corp. 

22  Meadow  Drive,  Stony  Brook,  New  York  11790 

SUMMARY 

An  iterative  algorithm  for  solving  nonlinear  inverse  problems  in  remote 

sensing  of  ocean  density  profiles  by  using  acoustic  pulses  is  developed. 

The  basic  idea  of  this  new  algorithm  is  that  first,  the  original  pulse 

problem  in  the  time-domain  is  reduced  to  a continuous  wave  problem  in 

frequency-domain  and  then  the  nonlinear  inverse  problem  in  frequency-domain 

is  solved  by  a hybrid  of  a Newton-like  iterative  method.  Backus  and  Gilbert 

linear  inversion  technique,  and  the  finite  difference  method.  This  new 

computational  algorithm  is  tested  by  numerical  simulations  with  given  data 

from  ten  different  frequencies  and  is  found  to  give  excellent  results.  The 

effects  of  taking  data  from  various  frequency  ranges  and  of  the  contaminating 

instrument  and  ambient  noise  on  the  accuracy  and  efficiency  of  numerical 

computation  are  investigated.  It  is  found  that  the  low  frequency  data  are 

preferred  over  the  data  from  the  high  frequency  spectrum.  Under  favorable 

conditions,  the  maximum  pointwise  numerical  error  of  the  density  profile  p(x) 

is  less  than  5%  of  the  total  variation  of  the  density  profile, 

ptv  = iMax.  p(x)  - Min.  p(x)|  . Better  result  can  be  achieved  if  larger 
x x 

number  of  data  are  available  and  more  efforts  are  made  in  the  numerical 
computation.  The  velocity  profile  can  be  obtained  from  the  inferred  density 
profile  through  a known  equation  of  state. 


INTRODUCTION 


Ocean  density  profiles  can  be  computed  from  a finite  number  of  experi- 
mental data  obtained  through  acoustic  remote  sensing  techniques  as  opposed  to 
in  situ  measurement.  Often  this  type  of  remote  sensing  problems  can  be  form- 
ulated as  "improperly  posed"  linear  or  nonlinear  inverse  problems  of  partial 
differential  equations  in  mathematical  analysis,  and  usually  the  solution  of 
an  inverse  problem  is  not  unique  and  does  not  depend  continuously  on  the 
given  data.  To  solve  linear  or  nonlinear  inverse  problems  is  equivalent 
to  constructing  approximate  solutions  of  linear  or  nonlinear  operator  equa- 
tions from  inadequate  data  with  or  without  statistical  measurement  errors. 

While  the  computational  methods  for  solving  the  traditional  direct  problems 
of  partial  differential  equations,  e.g.,  the  finite  difference  method,  the 
finite  element  method,  their  hybrids,  etc.,  are  quite  well  developed,  the 
numerical  methods  for  solving  inverse  problems  ere  still  in  their  infancies. 
Hence  the  need  of  new  developments  in  the  computational  methods  for  solving 
inverse  problems  is  definitely  in  order. 

Recently,  several  iterative  algorithms  for  solving  inverse  problems  of 
nonlinear  operator  equations  in  Banach  spaces  have  been  developed  by  Chen  and 
Surmont  [l],[2],  and  they  have  been  used  for  solving  a nonlinear  radiative 
transfer  problem  in  the  remote  sensing  of  atmospheric  temperature  profiles  [3]. 
These  algorithms  are  basically  hybrids  of  Newton's  iterative  methods  in  Banach 
spaces  [4]  and  the  linear  inversion  technique  of  Backus  and  Gilbert  (B  S G) 

[5], [6],  [7].  Of  course,  many  other  linear  inversion  techniques,  e.  g.,  the 
regularization  method  of  Tihonov  [8],  the  Moore-Penrose  pscudoinverse  method 
[9],  etc.,  can  be  used  just  the  same.  But,  the  B & G method  is  ' 'erred 
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here  for  the  reason  that  it  not  only  provides  an  inversion  technique,-  but  also 
can  be  used  as  a diagnostic  tool  for  testing  the  intrinsic  resolution  of  a 
given  set  of  data  for  a given  problem  (or  mathematically,  provides  a sort  of 
pointwise  error  estimates).  More  recently,  Tsien  and  Chen  [lO]  has  intro- 
duced a hybrid  of  the  above-mentioned  iterative  algorithms  and  the  finite 
difference  method  to  infer  the  ocean  density  profiles  by  using  low  frequency 
acoustic  continuous  wave  (CW)  measurements  with  success;  the  frequency  range 
is  such  that  the  ocean  depth  is  in  the  range  of  one  to  four  wavelengths. 

Now,  a numerical  algorithm  for  solving  nonlinear  inverse  problem  in 
remote  sensing  of  density  profiles  of  a simple  ocean  model  by  using  acoustic 
pulses  is  developed  in  this  paper.  Here  the  adiabatic  sound  velocity  is 
assumed  to  be  proportional  to  the  inverse  square  root  of  the  density  function. 
The  basic  idea  of  this  algorithm  is  that  the  incident  and  reflected  acoustic 
pulses  are  assumed  to  be  known  at  the  surface  of  the  ocean;  through  Fourier  sine 
or  cosine  transform  with  respect  to  time  t,  the  original  pulse  problem  in 
time-domain  is  reduced  to  a CW  problem  in  frequency-domain.  Hence  the 
numerical  method  of  [lO]  can  be  applied.  This  new  computational  algorithm 
is  tested  by  numerical  simulation  and  it  is  found  to  give  excellent  results. 

The  effects  of  taking  data  from  various  ranges  of  the  frequency  spectrum 
as  the  input  to  the  numerical  algorithm  on  the  accuracy  and  efficiency  of 
computations  are  thoroughly  studied;  the  entire  frequency  range  is  such  that 
the  ocean  depth  is  in  the  range  of  1 - 320  wavelengths.  Moreover,  the 
effects  of  the  contaminating  instrument  and  ambient  noise  on  the  accuracy 
of  the  numerical  results  or  the  intrinsic  resolution  of  a given  set  of  mea- 
sured data  are  carefully  investigated.  Finally,  a comprehensive  discussion 
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of  the  numerical  results  and  their  implication  in  actual  implementing  this 
computational  algorithm  is  given. 

It  should  be  pointed  out  that  there  are  two  basic  differences  between  the 
new  numerical  algorithm  here  and  the  conventional  approach  ([ll]  - [is])  for 
solving  inverse  Sturm-Liouvi 1 le  problems.  First,  a precise  knowledge  of 
complete  sets  of  eigenvalues  of  a particular  Sturm-Liouville  problem  and  its 
associated  problems  is  needed  to  solve  the  inverse  problem  uniquely  in  the 
conventional  approach,  but  no  such  knowledge  is  needed  in  the  new  numerical 
algorithm  here.  As  a matter  of  fact,  here  the  chioce  of  the  discrete  values 
in  frequency-domain  {u^},  i = 1,2,..., I,  (not  the  square  roots  of  the  eigen- 
values) is  completely  arbitrary.  Hence  they  can  be  chosen  for  practical  and 
computational  conveniences.  Secondly,  the  concept  of  the  uniqueness  of 
solutions  which  is  a major  concern  in  the  conventional  approach  is  replaced 
here  by  a weaker  but  more  practical  one,  "the  closeness  of  different  numerical 
solutions  to  a solution."  In  the  new  numerical  algorithm  here,  each  iterative 
calculation  with  different  initial  iterate  will  converge  to  a different  but 
unique  numerical  solution  [2],  nevertheless  this  non-uniqueness  poses  no  real 
difficulty  for  practical  problems.  If  several  different  numerical  solutions 
are  all  close  enough  to  the  exact  solution,  then  in  practical  sense  any  one  of 
them  will  be  an  acceptable  good  solution.  Hence  in  the  new  numerical 
algorithm  here,  the  major  concern  is  the  closeness  of  different  numerical 
solutions  to  the  exact  solution.  In  general,  the  "closeness"  is  characteri :^ed 
qualitatively  by  the  "spread  of  the  averaging  kernel",  i.  e.,  the  smaller  the 
spread,  the  closer  of  the  numerical  solution  to  the  actual  solution. 
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This  is  because  the  "spread"  characterizes  the  "closeness"  of  the  numerical 
solutions  of  the  individual  interate  - linearized  problem  and  the  ne;-/  numerical 
algorithm  here  can  be  shown  to  be  a Newton- like  iterative  method  whoso  accuracy 


is  in  direct  proportion  to  the  numerical  accuracy  in  computing  the  individual 
iterate  [2].  The  "spread"is  a by-product  in  the  process  of  computing  the 
numerical  solution.  The  quantitative  information  on  the  closeness  can  be 
obtained  by  performing  simple  numerical  simulations. 


FORMULATION  OF  PROBLEM 

The  hydrodynamic  structure  of  the  ocean  is  assumed  to  be  a compressible, 
inviscous  and  stratified  fluid  layer  above  a rigid  plane  (Fig.  1).  The 
linearized  governing  equations  in  remote  sensing  of  ocean  density  profiles 
by  using  acoustic  pulses  are  [l9] 


3u/3t 

= - p ^ 3p/3x  , 

0 < X < H 

3p/3t 

2 

= - pc  3u/3x  , 

0 < t , 

with  initial 

condi tions 

u(x,0)  = 0 , 

p(x,0)  = 0 , 

and  boundary 

c^'  ..i tions 

u(H,t)  = 0 , 

p(0,t)  = f(t)  , 

(1) 

(2) 

(3) 


where  p is  the  pressure,  u is  the  particle  velocity,  p is  the  equilibrium 


density,  c is  the  adiabatic  sound  velocity  related  to  p through  the  known 
equation  of  state  c = (p  K)"'^  [20] , K is  the  compressibility  coefficient 
which  is  assumed  to  be  constant,  and  f(t)  is  a pulse  function  with  a given 
duration  which  corresponds  the  incident  pulse. 

The  mathematical  problem  in  remote  sensing  of  ocean  density  profiles 
is  to  determine  the  density  function  from  Eqs.  (1),{2),(3)  and  3p/ax| 

x=o 


r 


To  transform  the  above-mentioned  pulse  problem  in  the  time-domain  to  its 
corresponding  CW  problem  in  the  frequency-domain,  Fourier  transform  with 
respect  to  t is  needed.  In  particular,  due  to  the  homogeneous  initial 
conditions  Fourier  sine  transform  is  applied  to  Eqs.  (1),(2)  and  (3). 

Upon  eliminating  the  Fourier  sine  transform  of  u , Eqs.  (1),(2)  and  (3)  leads 

to  a Sturm-Liouville  type  of  ordinary  differential  equation  for  £(x,u)),  the 

Fourier  sine  transform  of  p, 

d/dx  (p"^d£/dx)  + £ = 0 , 

£(0,o))  = f(o))  , (4) 

d£(x,u)/dx|^^^  = 0 , 

where  _f(u)  is  the  Fourier  sine  transform  of  f(t)  and  u is  the  angular  frequency. 

Upon  normalizing  K to  unity  and  introducing  the  dimensionless  variables, 
x'  = x/H,  £'=  £/£(w),  p'  = p/p(0),  and  w'  = wH/c(0),  Eq.  (4)  becomes 

-1  2 

d/dx  (p  d£/dx)  -i  u £ = 0 , 

£(0,0))  = 1 , (5)  j. 

d£(x,o))/dx|^^j  = 0 , 

where  p(0)  and  c(0)  are  the  values  of  p(x)  and  c(x)  at  x = 0 and  they  are 
assumed  to  be  known. 

Here  the  primes  are  dropped  for  convenience. 

Now,  the  inverse  problem  is  to  determine  p(x)  from  (5)  and  the  known 
values  of  d£(0,uj)/dx  for  a set  of  {o)^>,  i = 1,  2,...,I.  Although  Eq.  (5)  is 
a linear  ordinary  differential  equation,  the  inverse  proble’m  is  nonlinear 


in  nature. 
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NUMERICAL  ALGORITHM 

In  general,  a nonlinear  inverse  problem  can  be  solved  by  using  the 
iterative  algorithms  of  [l]  - [3].  For  the  nonlinear  inverse  problem 
here,  a variant  of  the  above-mentioned  iterative  algorithms  i-s  used  and  it 
is  presented  in  the  following. 

Let  -Efi+l  ~ ^ '^£n  ’ ^n+1  ” ^n  ^ *^*^n  * ^ ~ 0,1,2,....  (6) 

be  the  (n+l)th  iterates  of  £ and  p such  that  1 5£p  I Ifinl’  |pnl’ 

and  6pp(0)  = 0.  The  £^  and  are  the  initial  guess.  Upon  substituting 

p 2 

(6)  into  (5)  and  neglecting  terms  of  0{(6£^)  },  0{(5pp)  } and  higher,  one 

obtains  „ 

d/dx(p'^d£^/dx)  + 0,  £^  = 0 , ■ 

£^(0,w)  = 1 , dp^(x,a))/dx|^^^  = 0 , 

and  „ P 

d/dx(p'^d5£^ydx)  + w 62^  = d/dx(p^  6p,.,  dp^^/dx)  , 

(8) 

5£^(0,io)  = 0 , d6£^(x,co)/dx|^^^  = 0 . 

By  using  the  method  of  Green's  function  [21],  after  considerable  amount 
of  manipulation  one  obtains  a Fredholm  integral  equation  of  the  first  kind 
which  relates  Pp(x)  to  dfip^(x,u)/dx |^^q. 


/q  {p~^(x')  dp^(x' ,to)/dx' 6pj^(x')  dx'  = d5p^(x,o))/dx|^^g.  (9) 

Now  the  inverse  problem  is  reduced  to  for  each  iteration  the  determination 
of  Spp(x)  from  a set  of  measured  data  {d5£(0,u;'^)/dxl , i - 1,2,3,...,!.  These 
measured  data  are  used  as  the  input  data  to  the  right-hand  side  of  Eg.  (9) 
for  the  purpose  of  accelerating  the  rate  of  convergence.  This  inverse 
problem  can  be  solved  by  the  linear  inversion  technique  of  B & G [b]  - [7]. 
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It  starts  with  6pj^(x)  = T.  a-j,-,(x)  d5g(0,o)^- )/dx 

1 = 1 


(10) 


Upon  substituting  Eq.  (10)  into  Eq.  (9),  one  obtains  the  integral  expression. 


1 

'5Pn(x)  = / A (x,  x')op  (x')  dx'  , 

0 


v/here 


I 


A^(x,  x’)  = a.^(x)  {p|j^(x’)  d£^(x‘ ,ui)/dx'}'^  , 


m2 


(11) 

(12) 


is  known  as  the  " 
condition, 

The  "spread  of 


averaging  kernel"  and  is  subjected  to  the  normalization 

1 


/ A (x,  X')  dx'  = 
0 " 

from  x"  is  defined  by 


1 

Q„(x,A  ) : 12  / (x-x')‘ 

0 


1 . 


A^(x,  X')  dx'  . 


(13) 


(14) 


If  the  measured  data  are  errorless,  the  set  of  unknov/n  functions  {a^^(x)l, 
i = 1,2,3,...,!,  should  be  chosen  such  that  the  averaging  kernel  A^  resembles 
the  Dirac  delta  function  5(x'  - x)  most  closely  or  equivalently,  the  positive- 
definite  quadratic  form  of  {a^-p(x)},  i = 1,2,3,...,!,  Qj^(x,A^)  is  minimized 
subject  to  the  constraint  (13).  Hence  the  method  of  Lagrange  multipliers 
is  used. 


If  the  measured  data  contain  errors,  then  the  variance  of  tiie  density 
2 

error  Op(x)  incurred  at  x due  to  random  measurement  errors  can  be  found 
from  Eq.  (10)  to  be 

o2(x)  = A„(x)-E-A^(x)  , (15) 

where  the  vector  A^^(x)  = {a^^(x)},  i = 1,2,3,...,!,  E is  the  covariance  tensor 
for  the  random  errors  of  {d£(0,w^ )/dx} , i = 1,2,3,...,!,  and  the  means  of  the 
random  errors  are  assumed  to  be  zero. 
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Ideally,  one  would  like  to  be  able  to  choose  A|j(x)  such  that  both  o^(x) 
and  Q^(x,Aj^)  are  minimized.  However,  this  cannot  be  done,  but  it  is  possible 
to  minimize  a linear  combination  of  Cp(x)  and  Q^(x,A^), 

= (1  - s)Qn(x,A,^)  + s o^(x)  . (16) 

By  varying  the  parameter  s between  zero  and  unity,  the  emphasis  can  be 
shifted  from  minimization  of  the  spread  to  minimization  of  the  error.  Thus, 
there  is  a tradeoff  between  intrinsic  resolution  and  the  accuracy  in  the 
presence  of  measurement  errors.  The  best  choice  for  s must  be  determined  by 
the  particular  problem. 

Now,  for  a given  value  of  s the  minimization  of  Gp  subject  to  the 
constraint  Eq.  (13)  gives 

Ap(x)  = W'^(x)-Bp  / Bp-W-l(x)-Bp  , (17) 

where  the  matrix  Wp(x)  is  defined  by 

'X/  'V  'V 

Wp(x)  = (1  - s)Qp(x)  + s E (18) 

■'V/ 

with  the  elements  of  Qp(x), 

qijn(x)  = 12  / (x-x')2{p^l(x')d|^(x',ojT)/dx'}- 

{Pp^(x' )d£p(x' .Mj)/dx' } dx',  (19) 

i,  j = 1,2,3,...  ,I 

and  the  components  of  Bp, 

1 ' 
b^p  = / {p‘^(x‘)  dp^,(x‘ ,u^)/dx’ } dx'  , i = 1,2,3,...,!.  (20) 

0 

Finally  for  each  iterate,  once  A^(x)  is  known  from  Eq.  (17),  60p(x)  can  be 
calculated  from  Eq.  (10)._ 
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PROCEDURE  FOR  NUMERICAL  SIMULATION 

In  order  to  test  the  feasibility  and  to  study  the  general  characteristics 
of  the  new  computational  algoritiim  without  the  real  measurement  data,  a 
numerical  simulation  must  be  carried  out.  The  procedure  is: 

I.  First,  a typical  density  profile  of  the  ocean  [20]  p*(x)  is  assumed 

to  be  the  correct  density  function.  The  corresponding  pressures 
{£*(x,oj^)}  are  obtained  from  Eq.  (5)  for  a set  of  angular  frequencies 
{w^},  i = 1,2,3,...,!,  by  solving  the  two-point  boundary  value  problem. 

Then  by  a finite  difference  approximation,  {d£*(0,o)j )/dx} , i = 1,2,3,...,!, 
are  obtained  as  the  equivalence  of  the  measured  data. 

II.  The  measurement  errors  can  be  simulated  by  multiplying  the  data, 
{d£*(0,w,;  )/dx} , i = 1,2,3,...,!,  by  corresponding  factors, {(1  + k R^)}, 
i = 1,2,3,...,!,  where  k is  the  amplitude  factor  (k  = 0.001  used  here) 
and  R^'s  are  the  random  numbers  generated  by  a random  number  generator 
sub-routine  [22]. 

III.  Next,  Pq(x)  is  assumed.  Then  {£q(x,w^-)1  and  {d£^(0,u-j )/dx} , t 
i = 1,2,3,...,!,  are  determined  from  Eq.  (7)  by  using  the  same  numerical 
methods  as  those  of  Step  I.  Next,  the  non-homogeneous  term  of  Eq.  (9) 
can  be  computed  as  df£^(0,a)^- )/dx  = (1  + kRj)d£*(0,wj)/dx  - d£^(0,uj-j  )/dx. 

IV.  Finally,  S^Pq(x)  is  obtained  from  Eq.  (9)  and  the  data 
{d6£^(0,u^ )/dx) , i = 1,2,3,...,!,  by  using  the  linear  inversion  techni- 
que of  B S G [5]  - [7],  described  in  the  previous  section.  Hence 
from  Eq.  (6),  p^(x)  is  obtained. 

Other  than  the  truncation,  round-off  and  numerical  integration  errors, 
the  L2  norm  ||p*(x)  - p^(x)||  can  be  used  as  a criterion  for  e^aluatitig  the 
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perfornance  of  the  numerical  algorithm  consisting  of  steps  III  and  IV. 
Although  only  continuous  functions  are  considered,  I2  norm  is  used  here 
because  of  its  integrated  effect.  If  the  desired  accuracy  is  not  met,  then 
one  can  repeat  steps  III  and  IV  until  the  desired  accuracy  is  achieved. 

The  proof  of  the  convergence  of  this  iteration  method  will  be  shown 
elsewhere. 

Steps  I and  III  involve  solutions  of  the  two-point  boundary  value  problem 
of  the  Sturm-Liouvil le  type  of  ordinary  differential  equations.  In  general, 
a two-point  boundary  value  problem  can  be  solved  by  using  shooting  methods 
[23]  with  the  help  of  a sophisticated  ordinary  differential  equations  solver, 
GEAR  algorithm  [24], [25].  However,  in  the  range  of  u = 1 to  w = 30,  the. 

GEAR  algorithm  is  not  as  efficient  as  the  simpler  finite  difference  method 
with  second  order  accuracy,  e.  g.,  the  center  difference  approximation. 

In  the  range  of  w ^1,090  where  the  oscillation  of  the  solution  is  extremely 
high,  again  the  GEAR  algorithm  becomes  inefficient.  It  is  easier  to  obtain 
the  asymptotic  solution  of  the  two-point  boundary  value  problem  by  using  the 
WKBJ  method  [26] , [27] (Appendix  A). 

NUMERICAL  RESULTS  WITH  ERRORLESS  DATA 

Numerical  simulations  are  carried  out  for  several  p (x)  with  errorless 

0 

data  in  the  range  of  angular  frequency,  5 £ co  <_  2,000.  In  performing 
Stop  IV,  Simpson's  rule  is  used  in  the  numerical  integration.  The  numerical 
results  and  their  relevant  information  are  given  in  Table  1 and  Figs.  2 - '4. 

r\ 

hree  typical  oscillatory  {.-'[j^(x)  dp^(x , )/dx  r are  plotted  in  Fig.  10. 

The  averaging  kernels  A(x^,  x)  at  four  different  x^^  are  shown  in  Fig.  11. 


i 

r The  primes  of  the  independent  and  dependent  variables  arc  reinstalled  in  all 


tables  and  figures. 


■.n 

Fig. 

P'o(x’) 

n,  No.  of 
Iteration 

Method  Used 
in  Steps  I 
& III 

1 

I 1 3, 

t 1 f 

1 

i 

2 

6.98,  8.72,  10.12,  11.86, 
13.26,  15.00,  16.40,  18.14, 
19.54,  21.28 

Smooth. 
= 15.175 

3 

Impl icit 
Fini te 
Difference 
Ax’=n. 00625 

2.1802  1 

1 

i 

1 

1 

1 

Same 

as  Above 

7.igzaq, 
l|p'*-Pol|2 
= 100.13 

4 

Same 

as  Above 

1 

8.9105  1 

I 

. I 

1 

Same 

as  Above 

Smooth. 
= 1,340.4 

Same 

as  Above 

1 

10.009  i 

i 

1 

5 

100,  no,  no,  i3o,  i4o,  iso, 
160,  170,  180,  190 

Same  as 
Above 

9 

Shooting 

Method, 

"GEAR" 

23.121  ' 

1 
! 

6 

210,  215,  220,  225,  230,  310, 
315,  320,  325,  330 

Same  as 
Above 

2 

Same 

as  Above 

25.660  ! 

i 

■ 

410,  415,  420,  425,  430,  435, 
440,  445,  450,  455 

Same  as 
Above 

2 

Same 

as  Above 

54.434 

8 

10,  20,  30,  40,  50,  460,  550, 
640,  730,  820 

Same  as 
Above 

2 

Same 

as  Above 

21.241  : 

1 

9 

1100,  1200,  1300,  1400,  1500, 
1600,  1700,  1800,  1900,  2000 

Same  as 
Above 

2 

WKBJ  Method 

1 

601.92  I 

j 

TABLE  1 . Squares  of  norms  are  expressed  in 
units  of  10'^. 
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NUMERICAL  RESULTS  WITH  ERRONEOUS  DATA 

In  this  section,  the  effects  of  the  conta:.iinatinq  instrument  and 
ambient  noise  on  the  accuracy  of  the  numerical  results  or  tiie  intrinsic 
resolution  of  a given  set  of  measured  data  are  studied.  The  measurement 
errors  are  simulated  by  the  procedure  described  in  Step  II.  The  co- 
variance  tensor  of  ti)e  measurement  errors  E is  assumed  to  be  a diagonal 
matrix  with  elements 

e..  = {0.001  R^  d£*(0,ui)/dx}^  , i = 1,2,3,. ...I  . (21) 

In  consistence  with  the  study  of  the  effects  of  errorless  data  from  various 
ranges  of  the  frequency  spectrum  on  the  accuracy  of  the  numerical  results  in 
the  previous  section,  the  example  of  Fig.  4 is  used  for  the  study  iiere. 

The  effect  of  the  weight  parameter  s on  the  accuracy  of  the  numerical 
result  is  shown  in  Fig.  12  and  the  numerical  result  P2(>'')  for  the  case  s = 0 
is  plotted  in  Fig.  13.  The  tradeoff  curves  (02  vs.  s curves)  for  many 
ocean  depth  x^  = 0.00625,  X2  = 0.10625,  x^  = 0.20625,  x^  = 0.30625, 
x^  = 0.40625,  X0  = 0.50525,  Xy  = 0.60625,  Xg  = 0.70625,  Xg  = 0.80625, 
and  XjQ  = 0.90625  are  shown  in  Fig.  14.  Finally,  the  averaging  kernels 
corresponding  to  s = 0 at  four  different  ocean  depth  are  plotted  in  Fig.  15. 
Again  the  primes  of  the  independent  and  dependent  variables  are  reinstalled 
in  all  figures. 


DISCUSSION 


For  the  case  of  errorless  data,  all  of  the  numerical  results  in  Table  1 
and  Figs.  2-9  indicate  that  cur  numerical  algorithm  has  better  resolution 
near  the  ocean  surface  than  that  near  the  ocean  bottoni,  i.  e.  the  numerical 
solutions  in  general  are  more  accurate  in  the  upper  region  of  the  ocean  than 
other  regions.  This  conclusion  can  also  be  drawn  by  carefully  examining  the 
averaging  kernels  in  Fig.  11  where  the  averaging  kernels  resemble  the  Dirac 
delta-function  most  near  the  ocean  surface. 

The  numerical  results  in  Table  1 and  Figs.  2-4  suggest  that  the  closci' 
the  initial  guess  Pq(x)  to  the  exact  solution  p*(x),  the  more  accurate  the 
numerical  solution  will  be.  However,  it  can  be  shown  that  due  to  the  small- 
ness of  the  number  (inadequate}  of  data,  the  limit  of  the  iterates  differs 
from  the  exact  solution  and  deperids  on  the  initial  iterate.  The  detail  of 
this  has  been  discussed  in  [2].  Hence  in  the  real  computation,  the  initial 
guess  should  be  made  as  close  to  the  exact  proMle  as  possible. 

The  numerical  results  in  Table  1 and  Figs.  4-9  exhibit  the  fact  that 
data  from  lov/er  frequency  ranges  give  better  numerical  accuracy  than  those 
from  higher  frequency  ranges.  This  outcome  is  not  surprising  at  all, 
because  in  the  steps  I and  III  of  our  numerical  simulations  the  ordinary 
differential  equations  solvers  used  here  are  inefficient  in  the  high 
frequency  range.  At  the  same  time,  in  the  step  IV  of  our  numerical 
simulations  the  kernel  of  the  Fredholm  integral  equation  of  the  first  kind 
Eq.  (9)  is  much  less  oscillatory  in  the  lower  frequency  range  than  that  in 
the  higher  frequency  range  (Fig.  10).  This  means  that  the  numerical  solutioti 

of  Eq.  (9)  obtained  by  using  the  same  size  of  the  integration  sub-interval 
Ax  and  the  same  quadrature  formula  is  more  accurate  in  the  lower  frequency 


range  than  in  the  higher  frequency  range;  large  reduction  of  the  size  of 
AX  for  the  case  of  high  frequency  will  not  only  make  the  quadrature  fonnula 
more  inefficient  but  also  will  greatly  increase  the  globe  roundoff  error 
in  the  numerical  computation.  Hence  in  actual  computation,  low  frequency 
data  are  much  more  preferred  than  those  in  high  frequency  range.  Moreover, 
there  are  several  other  advantages  for  using  data  from  the  low  frequency 
range  and  they  will  be  discussed  later. 

For  the  case  of  erroneous  data.  Fig.  12  shows  that  the  accuracy  of  the 
numerical  result  decreases  as  the  parameter  s increases.  Even  when  s = 0 
(Fig.  13)  corresponding  to  t!ie  minimization  of  the  spread  alone,  the  accuracy 
is  poorer  than  the  analogous  case  with  errorless  data.  This  moans  that 
no  matter  how  small  the  contaminating  instrument  and  ambient  noise  can  be, 
they  will  lower  the  performance  of  the  numerical  algorithm. 

The  gross  behavior  of  the  tradeoff  curves  (Fig.  14)  is  the  same  for 
all  X.  The  random  error  decreases  monotonically  with  increasing  spread. 

The  decrease  is  rapid  at  small  values  of  the  spread,  followed  by  a rather 
abrupt  leveling  off  with  increasing  spread.  The  end  point  of  the  curve 
corresponding  to  minimum  spread,  and  maximum  random  error  corresponds  to  s - 0 
while  the  opposite  end  of  the  curve  corresponds  to  s = 1.  In  actual  co;r.- 
putation,  a happy  compromising  choice  of  s is  in  the  range  of  s corresponding 
to  the  neighborhood  of  the  bottom  of  the  steep  portion  of  the  tradeoff 

curves  (0. 5<( 1-s ) • 10^;$! ) . In  this  way,  the  effect  of  random  error  is  almost 

minimized,  while  the  spread  is  also  near  its  minimum.  Moreover,  upon  examin- 
ing the  tradeoff  curves  for  different  x,  one  finds  that  in  contrary  to  the 

case  of  errorless  data  where  the  nunserical  algorithm  has  the  best  resolution 
near  the  top  of  the  ocean  layer,  the  random  errors  in  the  data  have  shifted 
the  maximum  resolution  of  the  nun.erical  algorithm  to  somewhere  in  the  middle 
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of  the  ocean  layer.  This  fact  has  also  been  born  out  from  the  study  of  the 
averaging  kernels  for  different  x in  Fig.  15. 

In  real  situation,  there  are  three  major  reasons  favoring  the  use 
of  the  low  frequency  data  other  than  the  accuracy  of  the  numerical  algorithm. 
First,  the  low  frequency  component  acoustic  waves  arc  less  sensitive  to  the 
ocean  bottom  irregularities  and  other  scattering  marine  objects  tlian  tfie 
high  frequency  component  waves.  For  example,  if  the  ocean  depth  is  5,003  n 
and  the  average  sound  velocity  is  1,500  rn/sec,  then  the  normalised  angular 
frequency  range  5 £ w'  _<  50  corresponds  to  the  real  frequency  range 
0.239  Hz  f £ 2.387  Hz  with  wavelength  X in  the  range  62S  m < X <_  6,28d  m . 
Hence  it  is  conservative  to  say  that  our  numerical  algorithm  with  measured 
data  in  the  above-mentioned  low  frequency  range  will  not  be  sensitive  to  any 
irregularity  with  size  less  than  500  m . 

Secondly,  the  low  frequency  component  acoustic  waves  experience  less 
damping  due  to  viscosity  when  propagate  through  a fixed  distance  in  the  ocean 
than  the  high  frequency  waves,  so  the  received  reflected  signals  will  not  be 
drowned  out  by  the  contaminating  instrument  and  ambient  noise. 

Thirdly,  low  frequency  components  of  both  incident  and  reflected  acoustic 
pulses  contain  more  energy  than  the  high  frequency  components.  Hence  in 
practice,  the  adverse  effects  of  the  contaminating  instrument  and  ambient 
noise  are  minimized.  To  justify  the  statement  about  the  spectral  energy, 
a numerical  simulation  of  a short  acoustic  pulse  propagating  through  an  ocean 

layer  with  the  exact  density  profile  o*(x)  is  performed  in  Appendix  B . 

This  third  advantage  of  low  over  high  frequencies  can  be  ovf>-come  by  proper 

choice  of  hardware  designs,  e.g.  transducers. 

To  obtain  a more  realistic  but  conservative  evaluation  of  the  per- 


i 
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evaluation  criterion.  Upon  introducing  the  total  variation  of  the  density 

profile,  ; |Max.  p(x)  - Min.  p(x)|,  which  equals  to  0.00385p(0)  here 
o<_x<H  O'^x^H 

and  is  less  than  0.004p(0)  in  any  real  situation  [20]  , our  numerical  results 


are  given  in  Table  2 . 


Fig. 

Max.  |p*(x)-  p (x) j 
o<x<H ' 

Max.  Ip ’'■(x)- 

o<x<H  " 

Data 

2 

0.13  pj.y 

0.05 

errorless 

3 

0.64  pt^ 

0.13  p^y 

errorless 

m 

1-00  Ptv 

0.16  p,. 

errorless 

13 

1.00  p^^ 

0.18 

with  random  errors 

Table  2. 


(The  numerical  results  indicate  that  if  the  low  frequency  data  are  errorless  3 

and  the  initial  density  profile  is  close  enough  to  the  true  profile,  the  ' 

maximum  possible  pointwise  numerical  error  is  less  than  5"  of  the  total  ■ 

variation  of  the  density  profile.  The  accuracy  of  the  numerical  algorithm  i 

can  be  greatly  improved  if  there  are  more  measured  data  available  and  im- 
provements in  the  numerical  accuracy  are  made  for  each  individual  step  of 
the  numerical  algorithm.  However,  in  this  way  the  expense  incurred  in  the 
numerical  computation  will  be  greatly  increased.  Hence  in  practice,  a -j 

priori  decision  on  how  to  tradeoff  must  be  made  for  each  [jroblem. 

I ; 
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APPfNDIX  A 


For  (0  >>  100,  the  solution  of  the  two-point  boundary  value  problem  Fq.  (7) 
becomes  highly  oscillatory.  Hence  the  GEAR  algorithm  becomes  inefficient. 

An  excellent  method  to  circumvent  this  difficulty  is  to  construct  the  asymptotic 
solution  of  tlic  two-point  boundary  value  problem  for  large  to  by  using  the  WKBd 
method  [26], [27].  Again  in  this  and  next  appendixes,  the  primes  are  dropped 
for  convenience. 

First,  upon  setting  v,^(x,oO  = n^(x,to)/p[f(x) , Eq.  (7)  is  transformed  into 
v[](x,to)  + Jn(x,to)  v^(x,u))  = 0 , 

v^(0,to)  = 1 , (Al) 

v'(l,ut)  + p-^(l)  p;,(l)  vJl,to)  = 0 , 

where  d^(x,w)  = ‘vPn^(x)p[](x)  - (3/4 ) (p"^(x)pj](x) (A2) 
and  " ' " ^ d/dx.  If  Pp(x)  f 0 for  all  x and  Pp(x)  is  continuous  which  is  the 
case  here,  the  the  asymptotic  solution  of  Eq.  (Al)  is  given  by 


Vp(x,w)  = J~'^(x,tj){ap{(o)  cos  <t^(x,tj)  + bp(to)  sin  4.p(x,to)} 

where  <fp(x,.o))  = j J^(r,w)  dc  , 

^ 0 
U 

ap(u)  = Jp(0,w)  , 


(A3) 

(A4) 

(A5) 


and  bp(to)  = { j)((l ,oj) Jh'(0,M)sin  fp(l,v)  + 1 ,„  )Jp '( 1 ,w)jJf(0,oOcos  fn(l,  ) 

- p-^(l)p,;(l)J-'^(l,w)j;f(0,uOcos  ::„(l,-.o)}  / (A6) 

{jJ,=(l„o)cos  f^d.to)  - d};(l,to)J-/(l,.;)  sin  f.„(l,o) 

+ '-^■^(l)r'(l)J‘'^(I,to)sin  t^d,..))  . 

Then  p^(x,w)  = v^(x,(o)i ',',(x)  is  the  asymptotic  solution  of  Fq.  (7)  and  it  is 
used  in  the  integrand  of  Eq.  (9). 
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APPENDIX  8 


Tho  propacjation  of  a short  acoustic  pulse  throuqh  the  ocean  layer  with  the 
density  profile  p*(x)  is  described  by  the  solution  of  the  following  initial 
boundary  value  problem  of  the  linear  hyj'orbolic  system, 


Du/at  = 

- p*'^(x)  3p/ax  , 

0 < X < 1 , 

(Bl) 

9p/Dt  = 

- 3u/jX  , 

0 < t , 

with  initial  conditions 

u(x,0)  = p(x,0)  = 0 

9 

(B2) 

and  boundary  conditions 

u(l,t)  = 0 , 

(B3) 

p(0,t)  = 

100  n - cos(200^rt)} 

, 0 < t < 0.01 

= 

0 , 

t > 0.01  . 

To  damp  put  the  stable  numerical  oscillations  (Gibb's  phenomenon)  in  the 

neighborhood  where  the  solution  has  either  a discontinuity  or  a steep  gradient, 

2 2 

an  artificial  viscosity  term  va  u/sx*^  is  added  to  the  right  hand  side  of  the 
first  equation  of  Eq.  (Bl)i28j.  Along  with  it,  an  extra  boundary  condition, 
3u/ax|j^_Q  = - 20,000-rr  sin(200mt)  for  0 < t £ 0.01  and  0 for  t > 0.01 
must  be  added  to  Eq.  (83). 

The  above  system  is  solved  numerically  by  using  a method  which  is  based 
on  the  method  of  lines  with  the  space  variable  x discretized  by  a finite  element 
collocation  scheme  and  the  integration  of  time  by  the  GEAR  algorithm  [29]. 

In  particular,  the  Fourier  cosine  transform  of  u(0,t)  is  computed  by  using  a 
fast  Fourier  transform  routine  [30]  and  is  plotted  in  Fig.  16.  It  is  obvious 
that  the  low  frequency  spectrum  contains  more  energy  than  the  high  frequency 
spectrum.  If  the  pulse  has  longer  duration,  then  this  phenomenon  will  be  even 


more  pronounced. 
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1 Ocean  layer  with  a rigid  bottom. 

2 Comparison  of  the  calculated  and  the  exact  normalised  density  profiles, 
with  smooth  p'(x'),  small  l|p'*  - pAlj,  and  errorless  data  in 

6.5  < w'  < 2i:5. 

3 Comparison  of  the  calculated  and  the  exact  normalised  density  profiles, 
with  zigzag  po(x'),  small i|p'*  - pA[j,  and  errorless  data  in 

6.5  < w'  < 21.5. 

4 Comparison  of  the  calculated  and  the  exact  normalized  density  profiles, 
with  smooth  po(x'),  large  i|p‘*  - Pol^^  ^nd  errorless  data  in 

6.5  < u)'  < 21.5. 

5 Comparison  of  the  calculated  and  the  exact  normalized  density  profiles, 

with  smooth  ,o(x'),  large  j|p'*  - Poll*  ^nd  errorless  data  in 
100<u'<190.  p'*;  Pi; P^ 

6 Comaprisen  of  the  calculated  and  the  exact  normalized  density  profiles, 

with  smooth  oi(x'),  large  ]jp'*  - pA||)  and  errorless  data  in 
210<u,'<  330. p'*;__ Pi; 

7 Comparison  of  the  calculated  and  the  exact  normalized  density  profiles, 
with  Sniooth  ^;i(x'),  large  i|p'*  - pill>  and  errorless  data  in 

410  < co'  < 455. p'*; pi; p^. 

8 Comparison  of  the  calculated  and  the  exact  normalized  density  profiles, 
with  smooth  pi(x'),  large  ||p'*  - piM>  and  errorless  data  in 

10  to'  ^ 820. p'*;  pi; p^. 

9 Comparison  of  the  calculated  and  the  exact  normalized  density  profiles, 
with  smooth  ri(x'),  largo  i|p'*  - piM»  and  errorless  data  in 

1100  < to'  £ 2000.  p'*; Pi; P^. 

10  Typical  {oi”^(x‘ ) dpj!,(x ' ,ui' )/dx ' )^  as  function  of  x‘  shown  for  three 
different  values  of  w'. 

11  Typical  averaging  kernel  A(;<i,  x')  as  function  of  x'  shown  for  four 
different  Xq  with  errorless  data  in  6.5  < a)'  < 21.5. 

12  With  erroneous  data  in  6.5  < w'  < 21.5.  the  numerical  error  of  the 

example  in  Fig.  4 , characterized  by  Mp'*  - shown  to  be  a 

monotonic  increasing  function  of  s. 

13  Comparison  of  tiie  calculated  and  tlie  exact  normalized  density  profiles 
of  the  example  in  Fig.  4 with  erroneous  data  in  6.5  < w'  •-  21.5 

and  s = 0. 
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Fig.  14  With  erroneous  data  in  6.5  < oi‘  < 21.5,  tradeoff  curves  1,  2,  3,..., 
and  10  correspondin']  to  Xq  = 0.00625,  0.10625,  0.20625,  0.30625, 
0.40525,  0.50625,  0.60625','  0.70625,  0.00625,  and  0.90G25  respectively 
shoi.'ti.  The  top  and  botto;:i  end  points,  correspond  to  s = 0 and  s = 1, 

respectively.  x denotes  s 1 - 10"''  and  © denotes 

s = 1 - 0.5x  10-7. 

Fig.  15  Typical  averaging  kernel  as  function  of  x'  shown  for  four  different 
Xq  with  erroneous  data  in  5.5  < w'  < 21.5  and  s = 0. 

Fig.  16  Fourier  cosine  transform  of  u(0,t)  plotted  as  function  of  u'. 
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